Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

mkdir -p ~/EvaluationM4M5/FASTQ 
mkdir -p ~/EvaluationM4M5/CLEANING
mkdir -p ~/EvaluationM4M5/MAPPING
mkdir -p ~/EvaluationM4M5/QC
cd ~/EvaluationM4M5
tree ~/EvaluationM4M5

Téléchargement des données brutes

Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]

module load  sra-tools


srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir FASTQ

Combien de reads sont présents dans les fichiers R1 et R2 ?

cd FASTQ

cat SRR10390685_1.fastq | echo $((`wc -l`/4))
cat SRR10390685_2.fastq | echo $((`wc -l`/4))

Les fichiers FASTQ contiennent 7066055(R1) + 7066055(R2) = 14132110 reads.

Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

cd ../MAPPING/

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz

Quelle est la taille de ce génome ?


zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | awk '/^>/{if (l!="") print l; print; l=0; next}{l+=length($0)}END{print l}' |paste - -

La taille de ce génome est de 4215606 paires de bases.

Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz

Combien de gènes sont connus pour ce génome ?

zgrep -c $'\tgene\t' GCF_000009045.1_ASM904v1_genomic.gff.gz

4448 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit

cd ../
module load fastqc


gzip FASTQ/SRR10390685_1.fastq
gzip FASTQ/SRR10390685_2.fastq

srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8

La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

  • Oui
  • Non

car R2 comme le montre a plus de base dans la zone rouge que R1

Lien vers le rapport MulitQC

Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

  • Oui
  • Non

car il était nécessaire d’enlever “adapter content”

Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?

La profondeur de séquençage est de : 2.119142 G (total bases before filtering fastp) / 4215606 (ref genome size) = 502.7 (rather good) X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.

module load fastp

srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json CLEANING/fastp.json

Les paramètres suivants ont été choisis :

Parametre Valeur Explication

Ces paramètres ont permis de conserver 13.554096M reads pairés, soit une perte de 95,90992427882319% des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].

cd MAPPING/
GCF_000009045.1_ASM904v1_genomic.fna.gz > GCF_000009045.1_ASM904v1_genomic.fasta

module load bwa
srun bwa index GCF_000009045.1_ASM904v1_genomic.fasta


srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fasta ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685_on_GCF_000009045.1_ASM904v1.sam



module load samtools

srun --cpus-per-task=8 samtools view --threads 8 SRR10390685_on_GCF_000009045.1_ASM904v1.sam -b > SRR10390685_on_GCF_000009045.1_ASM904v1.bam

srun samtools sort SRR10390685_on_GCF_000009045.1_ASM904v1.bam -o SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam

srun samtools index SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam

srun samtools idxstats SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.idxstats

srun samtools flagstat SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.flagstat


samtools faidx GCF_000009045.1_ASM904v1_genomic.fasta

samtools faidx CP031214.1.fasta

Combien de reads ne sont pas mappés ?

13.57 M - 12.83M = 0.74 M reads ne sont pas mappés. Data from MulitQC report

13571369 - 12826829 = 744540 reads ne sont pas mappés. Data from flagstat and idxstats

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:



zgrep trmNF  MAPPING/GCF_000009045.1_ASM904v1_genomic.gff.gz  | awk '$3=="gene"' > trmNF.gff3


module load bedtools

bedtools bamtobed -i MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam  >  MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.bed

bedtools intersect -a MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.bed -b trmNF.bed  -f 0.5   > result_trmNF_bed.txt
grep -c 'NC_000964.3'  result_trmNF_bed.txt

2797 reads chevauchent le gène d’intérêt.

MulitQC report


export LC_ALL=en_US.UTF-8
export LANG=en_US.UTF-8

module load multiqc
srun multiqc -d . -o .

Visualisation

Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

References

1. toolkit NS. NCBI SRA toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.
4. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997. 2013.
5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
7. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (IGV): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.
 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France